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Abstract 

We  consider  a  time-dependent  problem  of  scattering  by  an  obstacle  involving  the  solution 
of  the  two  dimensional  Maxwell’s  equations  in  the  exterior  of  a  domain  with  a  perfectly 
conducting  condition  on  the  boundary  of  this  domain.  We  propose  a  novel  fictitious  domain 
method  based  on  a  distributed  Lagrange  multiplier  technique  for  the  solution  of  this  problem. 
Perfectly  matched  layers  are  constructed  to  model  the  unbounded  problem.  Comparisons 
are  performed  with  the  finite  difference  scheme,  that  demonstrate  the  advantages  of  our 
fictitious  domain  method  over  the  staircase  approximation  of  the  finite  difference  method. 
We  conclude  that  our  distributed  multiplier  approach  is  a  simple,  effective  and  far  more 
accurate  alternative  to  the  popular  FDTD  method  for  solving  Maxwell’s  equations. 
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1  Introduction 


Electromagnetic  phenomena  play  an  important  role  in  modern  technology  in  different  areas 
such  as  advanced  mobile  information  systems,  the  design,  development,  integration,  and 
testing  of  antennas,  communication  signal  processing  and  many  more.  Applications  involve 
the  propagation  and  scattering  of  transient  electromagnetic  signals  such  as  in  aircraft  radar 
signature  analysis  or  the  nondestructive  testing  of  concrete  structures.  The  study  of  such 
applications  requires  the  ability  to  predict  different  kinds  of  electromagnetic  effects.  Some  of 
the  important  effects  include  the  radar  scattering  attributes  i.e.,  radar  cross  section  (RCS) 
of  different  complex  objects  such  as  airplanes  and  missiles,  the  propagation  of  pulses  through 
dispersive  media  such  as  soil  or  concrete  to  detect  pollutants  or  hidden  targets,  interaction 
of  electromagnetic  waves  with  biological  media,  the  interaction  of  antenna  elements  with 
aircrafts  and  ships,  and  many  more  [23]. 

The  complete  set  of  laws  for  time-varying  electromagnetic  phenomena  can  be  derived 
from  physical  concepts  such  as  electric  charges  and  current  density,  some  universal  laws, 
such  as  the  conservation  of  electric  charge,  Faraday’s  and  Ampere’s  laws,  and  constituent 
laws  which  are  characteristic  for  a  given  medium  [11].  These  laws  of  electromagnetism  are 
represented  by  Maxwell’s  equations  and  are  central  to  predictions  such  as  those  described 
in  the  paragraph  above.  There  are  many  different  techniques  available  for  solving  the  time- 
dependent  problem  of  scattering  by  an  obstacle,  including  finite  difference  and  finite  element 
methods.  In  [5]  we  introduced  a  fictitious  domain  method,  based  on  a  distributed  Lagrange 
multiplier,  for  the  solution  of  the  two-dimensional  scalar  wave  equation  with  a  Dirichlet 
condition  on  the  boundary  of  an  obstacle.  In  this  paper  we  will  discuss  how  the  distributed 
Lagrange  multiplier  fictitious  domain  formulation  employed  in  [5]  can  be  applied  to  the  case 
of  the  two-dimensional  TM  mode  of  Maxwell’s  equations. 

A  fictitious  domain  method  is  a  technique  in  which  the  solution  to  a  given  problem  is 
obtained  by  extending  the  given  data  to  a  larger  but  simpler  shaped  domain,  containing  the 
original  domain,  and  solving  corresponding  equations  in  this  larger  fictitious  domain.  Let 
Cl  C  (d  =  1,  2,  3)  be  a  domain  which  contains  an  inclusion  u  as  shown  in  Figure  1.  We 
consider  solving  for  $  in  a  boundary  value  problem  of  the  type 

A(d>)  =  /,  in  Q\uj, 

Br($)  =  g0,  on  T  =  dfl,  (1.1) 

Beu($)  =  9i,  on  du, 

where  the  functions  /,  go,  <h  and  the  operators  A,  Br,  Bi)uJ ,  are  known.  In  a  fictitious  domain 
approach  we  replace  (1.1)  by  another  problem  of  the  type 


2 


n 


r 


Figure  1:  The  obstacle  u  embedded  inside  the  larger  domain  0. 


Find  0  defined  over  11  and  M1  a  measure  supported  by  du,  so  that 

(i)  A(<f>)  =  /  +  M7,  in  n, 

(ii)  Br((f>)  =  g0,  on  T  =  <911  (1.2) 

(iii)  £aw(0ln\fi>)  =  0i>  on  <9n;, 

where  the  operator  A  is  of  the  same  type  as  A  and  coincides  (in  some  sense)  with  donfi\w, 
/  is  some  extension  of  /  over  11  and  Bp,  and  BiJuJ  are  extensions  of  By,  and  Bg^,  respectively. 
If  we  choose  the  measure  M1  so  that  the  solution  of  (1.2,i,ii)  satisfies  relation  (1.2, iii)  then 
we  can  expect  to  have  0 |q\^  =  $.  If  H  has  a  simple  shape  such  as  a  rectangle  or  a  circle, 
for  example,  as  shown  in  Figure  1,  then  we  can  take  advantage  of  this  by  allowing  the  use 
of  uniform  finite  difference  grids  or  finite  element  meshes  and  hence  of  fast  solvers  for  the 
numerical  solution  of  the  finite  dimensional  systems  approximating  (1.1)  on  such  grids. 

Fictitious  domain  methods  can  be  traced  back  to  the  1960’s  to  Saulev  [26].  The  fictitious 
domain  method,  which  was  developed  to  handle  problems  with  complex  geometries  in  the 
stationary  case  [2,  16],  has  been  recently  applied  to  the  case  of  the  time  dependent  Maxwell’s 
equations  [7,  8,  9,  12].  In  all  these  cases,  a  boundary  Lagrange  multiplier  has  been  used 
to  enforce  the  Dirichlet  condition  on  the  boundary  of  the  obstacle.  In  [5]  we  compared  our 
distributed  multiplier  approach  to  the  boundary  multiplier  approach  for  the  one- dimensional 
scalar  wave  equation  and  observed  some  advantages  of  our  distributed  multiplier  formulation. 
We  refer  the  reader  to  [5]  for  more  details. 

The  advantage  of  our  distributed  multiplier  method  is  that  the  problem  in  the  fictitious 
domain  can  be  discretized  on  a  uniform  mesh,  independent  of  the  boundary  of  the  original 
domain,  thus  avoiding  the  time  consuming  construction  of  a  boundary  fitted  mesh  as  in  the 
finite  element  method.  (However,  there  are  some  classes  of  fictitious  domain  methods  that 
use  boundary  fitted  meshes  to  improve  accuracy  [19].)  At  the  same  time,  such  an  approach 
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is  more  accurate  than  the  staircase  approximation  of  the  finite  difference  method.  Thus,  the 
distributed  multiplier  approach  to  Maxwell’s  equations  presented  in  this  paper  provides  a 
simple  and  much  more  accurate  alternative  to  the  popular  FDTD  method  for  the  solution  of 
Maxwell’s  equations.  As  will  be  shown  here  the  approximation  of  the  distributed  multiplier 
approach  is  far  superior  to  the  staircase  approximation  of  the  FDTD  method.  In  addition  we 
also  demonstrate  the  ease  of  implementing  our  fictitious  domain  method  and  provide  details 
of  the  overhead  costs  which  are  minimal. 

We  consider  a  time-dependent  problem  of  scattering  by  the  obstacle  u.  We  are  interested 
in  calculating  the  field  in  the  exterior  of  u  C  with  d  =  2  or  d  =  3.  We  do  so  by  considering 
Maxwell’s  equations  in  free  space  with  a  Dirichlet  boundary  condition  on  du,  also  known 
as  a  perfectly  conducting  condition  (PEC).  The  evolution  problem  is  to  find  E  and  H  such 


that 


/ 

<9H 

(i) 

dE 

(ii) 

60  dT 

(iii) 

E  x  n 

,  (iv) 

E(x,  t 

in  \  u  x  (0,  T), 
in  Rd  \  u  x  (0,  T ), 


(1.3) 


3,  on  du  x  (0,  T ), 

))  =  E0(x),  and  H(x,  t  =  0)  =  H0(x),  in  M.d  \  u. 

In  (1.3),  E,  and  H  denote  the  electric  and  magnetic  fields,  respectively.  The  constants  e0,  and 
/To  denote  the  permittivity  and  permeability  of  free  space,  respectively.  The  speed  of  light  c 
is  given  to  be  c  =  1/y/eo ho.  Also,  n  is  the  unit  outward  normal  vector  to  the  boundary.  The 
initial  conditions  E0  and  H0  are  known  and  assumed  to  be  sufficiently  smooth.  This  is  an 
unbounded  problem.  One  of  the  ways  of  simulating  the  scattering  problem  in  an  unbounded 
domain  is  to  impose  an  absorbing  boundary  condition  on  the  boundary  of  the  truncated 
domain  0  enclosing  the  obstacle  u.  Hence,  we  consider  a  finite  domain  0,  and  we  impose  a 
first  order  (Silver-Muller)  absorbing  boundary  condition  on  the  (artificial)  boundary  F  =  <90. 
Thus,  our  time  dependent  problem  is  now  stated  as  an  evolution  problem  to  find  E  and  H 
such  that 


(i) 

(ii) 

(hi) 


dH  „  „ 

ho  +  VxE  —  0, 

dE 

eo—-VxH  =  0, 
dt 

E  x  n  =  0, 


m 


O  \u>  x  (0,  T), 


in  O  \  u  x  (0,  T), 
on  du  x  (0,  T), 


(1.4) 


(iv) 

l  (v) 


e0 


H  x  n  =  ,  / — n  x  (E  x  n), 

V  ho 

E(x,  t  =  0)  =  E0(x),  and  H(x,  t  =  0)  =  H0(x), 


on  T  x  (0,  T), 

in  O  \  u. 

We  will  also  consider  a  more  accurate  technique  called  a  perfectly  matched  layer  method  to 
simulate  such  unbounded  wave  propagation  problems  in  Section  7.  The  first  order  Silver- 
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Muller  boundary  conditions  on  T  model  the  electromagnetic  interactions  between  the  domain 
Q  and  the  exterior.  They  approximate  the  boundary  T  by  its  tangent  plane.  The  outgoing 
electromagnetic  plane  waves  which  propagate  normally  to  the  boundary  T  of  the  domain  0 
can  leave  freely  without  being  reflected  at  the  boundary  and  are  absorbed  at  the  boundary. 
The  Silver-Miiller  condition  on  T  x  (0,  T )  is  equivalent  to  the  Sommerfeld  radiation  held 
condition  for  the  Cartesian  held  components.  Applying  the  Silver-Miiller  conditions  at  a 
hnite  distance  from  a  spherical  scatterer  results  in  an  approximate  absorbing  boundary 
condition  which  is  exact  for  outgoing  spherical  waves  [22] . 

An  outline  of  the  paper  is  as  follows.  In  Section  2  we  present  a  distributed  Lagrange 
multiplier  formulation  for  the  problem  (1.4).  In  Section  3  we  present  wellposedness  results 
for  our  hctitious  domain  formulation  using  energy  methods.  In  Section  4  we  present  a  mixed 
hnite  element  formulation  for  the  (spatial)  discretization  of  the  hctitious  domain  problem 
and  we  present  an  iterative  method  for  its  solution  in  Section  5.  Decay  of  discrete  energies 
to  obtain  stability  results  are  proved  in  Section  6.  In  Section  7  we  replace  the  hrst  order 
absorbing  boundary  condition  by  perfectly  matched  layers  surrounding  the  computational 
domain.  Numerical  experiments  are  presented  to  validate  our  methods  in  Section  8  and  we 
conclude  with  an  outline  of  future  work  in  Section  9. 

2  A  Fictitious  Domain  Method  for  Maxwell’s  Equa¬ 
tions 

In  R3,  let  x  =  ( x,y,z ).  We  assume  that  neither  the  electromagnetic  held  excitation  nor 
the  modeled  geometry  has  any  variation  in  the  ^-direction.  That  is,  we  assume  that  all 
partial  derivatives  of  the  helds  with  respect  to  z  equal  zero,  and  the  structure  being  modeled 
extends  to  inhnity  in  the  z  direction  with  no  change  in  the  shape  or  position  of  its  transverse 
cross  section.  In  this  case  the  six  curl  equations  (1.4,i,ii)  can  be  decoupled  into  two  sets  of 
equations  each  involving  three  electromagnetic  held  vectors.  Let  E  =  (Ex,  Ey,  EZ)T ,  and  H  = 
(Hx,  Hy,  HZ)T  be  the  components  of  the  electric  and  magnetic  held  vectors,  respectively,  in  a 
Cartesian  coordinate  system.  In  the  TE  mode  the  electromagnetic  held  has  three  components 
Hz ,  Ex  and  Ey.  In  the  TM  mode  the  electromagnetic  held  has  the  three  components  Ez, 
Hx  and  Hy.  The  TE  and  TM  modes  are  decoupled  since  they  do  not  contain  any  common 
held  vector  components.  These  two  modes  are  completely  independent  for  structures  that  are 
composed  of  isotropic  materials  or  anisotropic  materials  in  which  the  off  diagonal  components 
in  the  constitutive  tensors  are  absent  [27]. 

We  will  consider  the  two-dimensional  TM  mode  of  Maxwell’s  equations.  Let  Q  now  be 
a  bounded  domain  of  R2,  with  x  =  ( x,y)T .  Let  H  =  (Hx,Hy)T  and  let  E  =  Ez.  Let 
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n  =  (nx,ny)T  be  the  unit  normal  vector,  and  let  us  define  the  unit  vector  t  =  (ny,—nx)J 
pointing  in  the  tangential  direction.  Then  system  (1.4)  becomes 

^ 

(i)  Ho— — h  curl  E  =  0,  in  fl  \  Co  x  (0,  T), 

ot 


dE 

(ii)  €q  — - curlH  =  0,  in  \  u  x  (0,  T), 

kJ  l 

(iii)  E  =  0,  on  du>  x  (0,  T), 

(iv)  H  •  t  =  on  T  x  (0,  T), 

V  E o 

(v)  E{x,  t  —  0)  =  Eq(x.),  and  H(x,  t  —  0)  =  H0(x),  in  \  Co. 


(2.1) 


In  the  above,  the  operator  denoted  by  curl,  is  a  linear  differential  operator,  which  is  defined 
as 

- »  fin  rin 

(2.2) 


^ v“ep'(a)' 


where,  V'(Vt)  is  the  space  of  distributions  on  O.  Similarly,  the  linear  differential  operator 
denoted  by  curl  is  defined  as 


curlv=^-^  V  v  =(vx,vv)e  V(a)2. 


(2.3) 


The  operator  curl  appears  as  the  (formal)  transpose  of  the  operator  curl  [10],  i.e., 

(curlv,  0)  =  (v,  curl0),  Vv  G  T>'(Q)2,  <f>  G  (2.4) 


In  the  case  of  the  two  dimensional  TM  mode,  the  PEC  condition  E  x  n  =  0,  on  du>  x  (0  ,T) 
translates  to 

E  =  EZ  =  0,  on  dui  x  (0,  T).  (2.5) 

We  assume  that  the  fields  (E,  H)  are  sufficiently  differentiable  in  time.  We  note  that  (2.1,iv) 
is  the  Silver-Miiller  condition  (1.4,iv)  for  the  TM  mode.  The  cross  product  H  x  n  can  be 
written  as  H  •  tz,  where  5  is  a  unit  vector  in  the  z  direction. 

We  employ  the  fictitious  domain  method  introduced  in  [5]  to  enforce  the  Dirichlet  bound¬ 
ary  condition  (2.1, iii)  on  the  boundary  du j  of  the  obstacle  uj.  The  Silver-Miiller  boundary 
condition  is  naturally  incorporated  into  the  weak  formulation  that  we  construct  by  integrat¬ 
ing  (by  parts)  the  equations  (2. 1  ,i,  ii)  over  the  domain  Q,  in  (2.1,ii).  Thus,  the  Silver-Miiller 
boundary  condition  does  not  have  to  be  enforced  in  the  functional  spaces.  Using  a  dis¬ 
tributed  Lagrange  multiplier  approach  problem  (2.1)  is  equivalent,  at  least  formally,  to  the 
variational  problem 


6 


Find  {E(-,t),  H(-,  t),  A(-,  t)}  £  H1^)  x  [L2(fl)]2  x  L2(u>)  such  that 


d 


(i)  no  —  /  H  ■  dx  +  /  curli?  •  \F  dx  =  0,  V  \F  e  [L2(fl)]z 
dt  Jq  Jn 


(ii)  eo-j-  f  E<p  dx  —  /'H-curl0dx 
dt  Jo  Jo 


£o 
ho  dr 


Ecf)  dr, 


+  /  A0d(u  =  O,  V  (f>£H\n), 

J  uj 

(iii)  j  Et  den  =  0,  V  r  G  L2{uj), 

Joj 

\  (iv)  E{x,t  —  0)  —  E0{x),  and  H(x,  t  —  0)  =  H0(x),  in  id, 
in  the  sense  that 


E  = 


E  on  Q,\u>, 
0  on  duj. 


;  h 


H  on  £l\u, 
0  on  duj. 


(2.6) 


(2.7) 


The  function  Eq  is  chosen  to  be  an  Hl  -  extension  of  Eq,  and  Ho  to  be  at  least  an  L2-extension 
of  H0.  Thus,  we  have 


£(x,  t  =  0)  = 


E'o(x)  on  Q\u, 


0 


on  uj. 


,  H(x,  t  =  0)  = 


H0  (x)  on  Q\u, 


0 


on  uj. 


(2.8) 


_ ^  QJJJ  QJJJ 

We  note  that,  for  E  £  L2(fl),  curlF1  =  (— — ,  —  — — )  £  \L2(fl) l2,  implies  that  both  the  partial 

oy  ox 

derivatives  of  E  must  be  in  L2(Q).  Hence  we  must  have  E  £  Hl(Vt).  In  succeeding  sections 
we  will,  however,  drop  the  ~  symbol  on  the  fields  E  and  H.  Thus,  the  system  (2.6)  will  read: 
Find  {E(-,t),H(-,t),\(-,t)}  £  H\Q)  x  [L2(Q)]2  x  L2(uj)  such  that 

/  (i)  hot  [  n-V  dx+  [cmlE-V  dx  =  0,  V  V  £  [L2(Sd)]2, 
dt  Jn  Jn 


d 


(ii)  e0—  Ecf)  dx-  /  H  •  curie/)  dx  +  J  —  /  E(\>  dr, 


dt 


IQ 


Mo  Jr 


+  /  \<j)  den  =  0,  V  (f)£H\Q), 

J  UJ 

(iii)  J  Et  den  =  0,  V  r  £  L2(uj), 

[  (iv)  F?(x,  t  —  0)  =  Eq(x.),  and  H(x,  t  =  0)  =  H0(x)  in  Q. 


(2.9) 


The  idea  behind  the  distributed  fictitious  domain  method  is  to  extend  the  electromagnetic 
solution  inside  the  obstacle  uj,  and  solve  Maxwell’s  equations  in  the  entire  domain  Cl.  The 
Dirichlet  condition  on  duj  is  enforced  via  the  introduction  of  a  Lagrange  multiplier  on  the 
entire  domain  uj.  In  [6,  12]  a  boundary  multiplier  fictitious  domain  method  is  introduced  for 
the  wave  equation,  and  for  Maxwell’s  equations. 
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3  Conservation  of  Energy 


In  this  section  we  derive  an  energy  identity  from  the  variational  formulation  (2.9).  The  energy 
identity  presented  below  guarantees  the  wellposedness  of  the  problem,  and  the  stability  of 
the  solution.  Let  (•,  •)  denote  the  L 2  inner  product  in  fl,  (•,  the  L2  inner  product  in  u>, 
and  (•,  -)r  the  L2  inner  product  on  T.  Also,  let  ||  •  ||?||  •  ||^?||  *  ||r  denote  the  corresponding 
norms. 


Theorem  1  The  system  (2.9)  verifies  the  energy  identity 


where  the  energy  £  is  defined  as 


1 


with 


T  =  -{e0||£||2  +  /r0||H||2}, 


\  1/2 

I2  dr) 


Thus ,  (3.1)  implies  that  the  energy  does  not  grow  over  time,  i.e., 

£(t)  <  £(0),  Vt  >  0. 

Proof.  Let  us  take  <f>  =  E  in  (2.9,ii).  We  obtain 

e0-y-  [  \E\2  dx  —  [  K-cmlE  dx  +  .A  [  |£|2dT  +  [  XE  den  =  0. 
dt  Jq  Jq  y  p0  J r  Ju 

Next,  we  take  SI/  =  H  in  (2.9,i).  With  this  choice  we  get 

Po  ——  [  |H|2dx+  [  curl E  •  H  dx  =  0. 

dt  Jn  Jn 

Adding  equations  (3.5)  and  (3.6)  we  have 

eo~r,  [  \E\2  +  pA  [  |H|2  dx  +  [  \E\ 2  dr+  [  XEdoj  =  0, 

at  Jci  at  Jq  y  p0  Jr  Ju 

which  can  be  rewritten  as 

~  (eo  ||£||2  +  p0  ||H||2)  +  \A  f  \E\ 2  dr+  [  XEduj  =  0. 

^  dt  V  Mo  Jr  Juj 

Taking  r  =  A  in  (2.9,iii)  we  obtain 


/ 

J  i 0 


EX  d uj  =  0. 


(3.1) 

(3.2) 

(3.3) 

(3.4) 

(3.5) 

(3.6) 

(3.7) 

(3.8) 

(3.9) 


Substituting  (3.9)  in  (3.8),  and  using  the  definition  of  the  energy  (3.2)  we  obtain  (3.1).  ■ 
Equation  (3.1)  implies  that  there  is  no  dissipation  of  the  waves  in  the  domain  £1.  This 
is  the  principle  of  conservation  of  energy  for  the  variational  formulation  (2.9)  for  Maxwell’s 
equation. 


Figure  2:  (left)  A  staircase  approximation  to  a  scattering  disk.  The  disk 
is  approximated  by  the  highlighted  nodal  points,  (right)  The  degrees 
of  freedom,  for  the  Lagrange  multiplier  A  in  the  fictitious  domain 
method,  in  the  case  of  a  scattering  disk.  The  mesh  ratio,  i.e.,  the  ratio  of 
the  step  size  chosen  on  the  obstacle  to  the  mesh  step  size,  is  about  1.3. 

4  Numerical  Discretization:  A  Mixed  Finite  Element 
Method 

A  very  popular  technique  used  to  solve  Maxwell’s  equations  is  the  finite  difference  time 
domain  method  (FDTD)  which  uses  a  rectangular  grid  and  an  explicit  scheme  in  time.  The 
degrees  of  freedom  of  the  electric  and  magnetic  field  are  staggered  in  space  and  time.  This 
method  is  computationally  very  efficient,  however  the  staircase  approximation  to  the  obstacle 
is  inaccurate,  and  it  leads  to  excessive  numerical  diffraction  when  the  obstacle  boundary  does 
not  fit  the  mesh,  as  seen  in  Figure  2  (left).  In  this  figure  the  scattering  obstacle  is  a  disk, 
and  is  approximated  by  the  darkened  nodal  points.  We  now  present  details  of  the  numerical 
approximation  of  problem  (2.9) 

4.1  Space  Discretization 

Let  now  be  a  union  of  rectangles  such  that  we  can  consider  a  regular  mesh  (%,)  with 
square  elements  (K)  of  edge  length  h  >  0  as  in  Figure  3.  The  approximation  space  for  the 
electric  field  E  is  chosen  to  be 

Uh  =  {(i>heH1{n) I  VAT  e  %,(j)h\K  e  Qi},  (4.1) 

where,  the  space  Q\  =  P\  \ .  The  basis  functions  for  E  have  unity  value  at  one  node  and  are 
zero  at  all  other  nodes.  Figure  3  shows  the  locations  for  the  degrees  of  freedom  for  both 
approximation  spaces. 
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Edge  2 


Figure  3:  A  sample  domain  element  K.  The  degrees  of  freedom  for  the 
electric  and  magnetic  field  are  staggered  in  space.  The  degrees  of  freedom 
for  the  electric  field  E  are  at  the  nodes  of  the  square.  The  degrees  of 
freedom  for  Hx  and  Hy  are  the  midpoints  of  edges  parallel  to  the  a;- axis 
and  y-axis,  respectively. 


For  the  approximation  of  the  magnetic  field  H  we  need  to  consider  a  space  Vh  that 
satisfies  the  property 

VxUhcVh  (4.2) 

in  order  to  ensure  a  well  posed  formulation,  see  [20,  21,  17]  for  sufficient  conditions  for 
convergence  of  mixed  methods.  To  this  end  we  choose 

Vh  =  g  [L2Q)]2\  VK  g  Th,  *h\K  e  RT[0]},  (4.3) 

where,  RT[q]  =  P±o  x  P0i>  is  the  lowest  order  Raviart  Thomas  space  [24]  and  for  k \ .  k2  G 
N  U  {0}, 

Pk^  =  {p(xi,x2)\p(xl,x2)  =  ^2  a^x \x32}. 

0<i<&i  0<j<k2 

The  basis  functions  for  Hx  and  Hy  have  unity  value  along  one  ey  or  ex  edge,  respectively, 
and  zero  over  all  other  edges  (see  Figure  3). 

Let  the  set  of  mesh  points  on  fi  be  defined  as 

E/i  =  {P  |  P  G  Cl,  P  is  a  vertex  of  Th}.  (4-4) 

Next,  we  define  the  set 

E h  =  {P\P  G  7),  d(P,  dcu)  >  h}  U  Discrete  set  of  points  belonging  to  dco.  (4.5) 
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The  points  on  dcu  are  typically  chosen  so  that  their  distance  is  of  the  order  of  h.  Using  the 
sets  defined  above,  we  now  define  the  set  A/,  of  the  Lagrange  multipliers  by 


Aft  —  {//ft  |  /ift  —  ^  /GpXp,  Up  £  (4-6) 

Pesg 

with  xp  the  characteristic  function  of  the  elementary  square  of  center  P  and  edge  length 
h\  we  clearly  have  ///,  ( P)  —  /ij>.  We  approximate  the  integrals  involving  the  distributed 
multiplier  by 


/ift  vh  d x 


h2  yUft(P)  vh(P),  Vnft  G  Vft,  V/ift  G  Aft. 

PGSg 


(4.7) 


Figure  2  illustrates  a  choice  for  the  set  in  the  case  of  a  scattering  disk.  The  ratio  of  the 
distance  between  points  on  the  circle,  denoted  by  hg^,  to  the  mesh  step  size,  h,  is  about  f.3. 
We  will  call  this  ratio  as  the  mesh  ratio.  In  numerical  experiments,  good  results  are  observed 
when  the  mesh  ratio  is  approximately  1.5  or  greater  [14].  Based  on  these  approximation 
spaces,  the  space  discrete  scheme  is  defined  as: 

Find  {Pft(-,t),Hft(-,f),  Aft(-,f)}  G  Uh  x  Vft  x  Ah  such  that 


(i)  H o-j:  [  Hft  •  \Fft  dx  +  [  curl  Eh  ■  ^fh  dx  =  0,  V^G  Vft, 

dt  Ja  Jn 

(h)  e0—  /  Eh<f>h  dx-  /  Hft  •  curl (ph  dx  +  J  —  /  Eh<j>h  dT, 

Jn  Jn 


(hi) 


dt  jq  j a 

+  /  Aft/^ft  da;  =  0,  V  ^ft  G  7/ft, 
J  UJ 

[  Ehrh  da;  =  0,  V  rh  G  Ah, 


Mo  Jr 


(4.8) 


{  (iv)  Eh(x,t  =  0)  =  P0,ft(x),  and  Hh(x,  t  =  0)  =  H0,ft(x)  in  ft. 


4.2  Time  Discretization 


For  the  time  discretization  we  use  a  centered  second  order  accurate  finite  difference  scheme 
in  which  the  electric  and  magnetic  field  components  are  staggered  in  time,  1/2  units  apart. 
For  example,  we  compute  the  magnetic  field  at  the  time  step  n  +  1/2  and  the  electric  field 
at  the  time  step  n  +  1.  On  the  interval  [0,  T],  let  At  =  T/N  be  the  time  step,  where  N  G  N. 
We  define  Vk  =  V/ft  =  kAt)  and  denote  tk  =  kAt,  for  k  G  Z.  Following  the  notation  in  [3] 
we  define  for  k  G  Z 

T/fc+1/2  _  TP&-1/2 

D*tVk  = - - - ,  (4.9) 


and 

—k  l/^+V2  +  vk~1/2 

vk  =  - - - .  (4.10) 
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Using  the  notations  and  definitions  developed  above,  the  fully  discrete  scheme  is  given  as 
For  n  —  0, 1, . . . ,  N  —  1,  on  the  interval  ( tn ,  tn+1),  given  E%,  H+1//2 
Find  {El+1,  H”+1/2,  A+1}  G  Uh  x  Vh  x  Ah  such  that 


(i) 

po(DAtKnh,*h)  +  (cm\El*h)  =  0:  V^eVfc, 

(ii) 

eo  (DAtEl+1'2^h) 

~  (H^+1/2,  curl  <ph)  +  ^ 

+(S"+1/2.^)r. 

1  Po 

+Ui)+1>  01+  —  0, 

V  4>h  G  W/j, 

(4.11) 

(iii) 

(B;+1,rh)„  =  0, 

n  ^  A./; , 

(iv) 

E«(x)  =  E0ih(x). 

and  H“1/2(x)  =  H0, hl 

x)  -  ^  curl£’o,/l(x)  in  D. 
Zp0 

We  will  use  quadrature  rules  for  the  calculation  of  all  integrals  involved  in  (4.8)  due  to  which 
we  obtain  diagonal  mass  matrices  and  an  explicit  scheme  in  time.  The  use  of  quadrature 
formulas  to  obtain  diagonal  mass  matrices  is  referred  to  as  mass-lumping.  Similarly,  we  use 
quadrature  rules  to  calculate  the  boundary  integrals  as  well.  In  this  case  system  (4.11)  in  the 
absence  of  the  distributed  Lagrange  multiplier,  reduces  to  the  FDTD  scheme  on  a  regular 
mesh. 


5  Iterative  Solution  of  the  Discrete  Problem 


Let  £+e|n+1  be  the  electric  field  solution  to  (4.11)  in  the  absence  of  the  distributed  multiplier 
(i.e.,  solution  to  the  FDTD  scheme).  Then  the  update  equations  for  the  scheme  (4.11)  can 
be  presented  as  follows.  For  an  interior  node  (l.  to)  we  have 


(i)  H, 


in+1/2 
xU,m+ 1/2 

n+1/2 


=  H 


in— 1/2 


x  U, m+1/2 
n- 1/2 


At 

Poll 

At 


/••A  tt  n-\-L/z  _  tt  n—i/A _ _/ 771711  77>7l|  \ 

W  ny  1+1/2, m  ~  ny  1+1/2, m  +  u  \l+l,m  &  \l,m b 


At 


Poh 

I  n+1/2 


-H, 


in+1/2 


(iii)  Iff  =  E\lm  +  ^Hy\Z lira  -  -.11-1/2 ,m 

°A  t 


) 


—AH, 


in+1/2 


_  tt  n+1/2  \ 

£lx  Lm- 1/2/’ 


e0h  y  X  lz»m+i/2  ±±x\l,m—l/2 


(5.1) 


where  for  a  solution  component  Vjf,  V | +  =  V/ (x  =  lh,y  =  mh ),  with  k,l,m  G  Z.  For  a 
node  on  the  boundary  F,  the  boundary  integral 


J  Ehn+1/2(f)hdr,  (5.2) 

will  contribute  terms  to  both  the  right  hand  side  and  the  left  hand  side  of  equation  (5.1, iii), 
as  this  term  involves  £++  which  is  unknown,  as  well  as  E%,  which  is  known.  In  this  case 
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(5.1,iii)  has  to  be  modified  as 


(iii)  ErCn  =  ”rn  +  ^SHlS1/2-  (5.3) 

In  the  above 

/ 1  ncAt\  _  ( 1  ncAt\ 

7  7  =  +  ah  J  ’ 

where,  for  an  interior  node  (5  =  1,  a  =  1,  k  =  0,  for  a  boundary  node  but  not  a  corner  node 
f3  =  2,  cr  =  2,  k  =  1,  and  for  a  boundary  corner  node  [3  =  4,  a  =  4,  k  =  1.  Also,  S'  is  the 
stiffness  matrix  associated  with  the  integral 


/  H™+1/2  cVl  <ph  dx,  G  Z4-  (5.5) 

Jn 

The  solution  E)}t  + 1  to  the  scheme  (4.11)  is  obtained  from  the  solution  £'^ee|n+1,  to  the  FDTD 
scheme  (including  absorbing  boundary  terms),  by  adjusting  for  the  Dirichlet  condition  on 
the  obstacle  uj  via  the  Lagrange  multiplier  Ajj+1.  Thus,  we  will  solve  a  system  of  the  form 
Find  (££+1,  \nh+1)  eUhxAh  such  that: 


'  DhEZ+1  +  spr1  =  £TI”+1. 
,  BkK+1  =  0. 


(5.6) 


where  Dh  is  the  lumped  mass  matrix  associated  to  the  integral  /  E^Ok  dx,  Vb/,.  G  Uh, 

Jn 

and  the  matrix  Bh  is  associated  with  with  the  integral  (4.7).  We  note  that  Dh  G  MiVxiV  is 
symmetric  positive  definite,  and  Bh  G  MMxJV(M  <<  N ).  We  use  the  Schur  Complement  of 
the  system  (5.6) 

(BkDpBTh)  A"+1  =  BhDpEC  |”+\  (5.7) 

to  solve  for  Ajj+1.  We  do  this  by  using  an  Uzawa-Type  conjugate  gradient  algorithm  in  the 
form  given  in  [15].  As  remarked  in  [5],  the  matrix  BhDjll  Bj  is  symmetric  and  positive 
definite,  a  property  that  is  related  to  a  uniform  discrete  inf-sup  condition  associated  with 
the  integral  (4.7).  We  refer  the  reader  to  [5,  14]  for  further  details  on  inf-sup  conditions  for 
mixed  problems. 


6  Decay  of  Discrete  Energies 

In  this  section  we  prove  a  discrete  energy  identity  based  on  the  discretized  fictitious  domain 
formulation  (4.11).  This  identity  indicates  the  stability  of  the  distributed  multiplier  formu¬ 
lation.  An  interesting  fact  to  note  is  that  the  stability  condition  (CFL)  is  the  same  as  in  the 
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case  of  the  problem  without  an  obstacle.  We  define  the  bilinear  form 

a(cf)h,  iph)  =  /  curl  •  curl  dx  ,  V  (<f>h,  iph )  eUhx  Uh, 
Jn 


and  the  operator  Kh'-lAh^  hf'h  by 

(fch<t>h,4’h)n  =  a(<l>h,4’h)- 

Let  /  be  the  identity  operator  on  Uh- 

Theorem  2  If  the  CFL  condition  cAt  <  h/V 2  is  satisfied,  then  the  operator 

c2  At2 

Sh  =  I-  —)Ch, 

defines  a  positive  quadratic  form,  the  expression 


{f‘ol|Hr1||2  +  £oK+1,5»i5+1)}, 


defines  a  discrete  energy,  and  system  (4.11)  verifies  the  energy  identity 

£%+1  =  ££-  At,[^\\Enh+1/2\\2r,  Vn  e  N,n  >  0.  (6.5) 

V  ho 

Thus,  (6.5)  implies  that  the  discrete  energy  does  not  grow  over  time,  i.e., 

£h+1  <£h,  Vn  >  0.  (6.6) 

Proof.  We  consider  the  mean  value  of  (4.1 1  ,i)  at  n  and  n  +  1,  with  SEu  =  jj«+1/2  ge^ 

Po(^hT+1/2,H^+1/2)  +  (^\Ehn+1/2,Unh+1/2)  =  0.  (6.7) 

Using  <fh  =  E^+1^2  in  (4.11, ii)  we  get 

£o(BA«K+1/2,B;+1/2)-(Hp1/2WlB;+1/2)+,/V|E+‘/2|lr+(Ar1,E;'+‘/2V  =  o.  (6.8) 

V  ho 

Taking  the  mean  value  of  (4. 1 1  ,iii)  at  n  and  n  +  1  and  by  taking  fit,  =  A)'+1.  we  have 

(eT1,2K+\  =  0-  (0.9) 


Adding  (6.7),  (6.8),  and  (6.9)  we  get 

MBIHT,+1/2,  Hr1/2)  +eo(CatBy1/2,£y+1/2)  =  -M\\E^+1/2\\l 


(6.10) 
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This  implies 


2A  t 


h: 


+3/2  j j  72+1/2 


)+€b||KflHS 


2A7 


H 


72+1/2  TT  72 — 1/2\ 


+ 


eolRII2}  - 


-ll£T+I/2||?.. 


Using  the  parallelogram  law  we  can  write 


Hn+3/2  jj-n+l^l  _  1  i 

h  )■*+  J  —  ^  I 


H 


■72+3  /  2 


+  H 


n+1/2  m2 


Similarly, 


TT^m/2  Tjfl  — 1/2 

+  )  *+ 


=  l|H, 

=  I|Ha 


=r= — ™+l  i|2 


|H 


72+1  1 1  2 

h  II* 


■72+3  /  2 


H 


n+1/2  i|  2 


n+1/2 


+  H 


n— 1/2  1 1 2 


il'H 


■n+1/2 


H 


n— 1/2  i|  2 


”"2  - +|£a,h;;||2. 


From  (6.12)  and  (4.1 1  ,i)  we  have 

1  f  /-w-T?i+ 3/2  tj-tz+1/2' 


H 


2  ^  Mo||H, 


1+1 1 12 


>  +  eol«+1||2} 

At2 


4/+o 


curlK+1H2  +  eolK+1|i 


=  4”+1- 


H 


i.”+1||2  +  €o(K+1.^K+1)} 


Similarly  we  can  show  that 

1 

2 


H 


72+1/2  tj-72— 1/2' 


+  eol|Bh"H2}=£; 


(6.11) 


(6.12) 


(6.13) 


(6.14) 


(6.15) 


Substituting  (6.14)  and  (6.15)  in  (6.11)  we  obtain  the  energy  identity  (6.5)  In  two  dimensions 
we  have, 

h2{JCh(f>h,(f>h)i?(n)  0 

suPfceu»^«ur<  ’ 

which  along  with  the  CFL  condition  implies  that 

c2At2 , 


(6.16) 


(++/+/i)l2(ST)  =  (4>h,  (7 - Z - K'h)<fih) L2 (Cl)  >  0,  V(j)h  eU h  \  {0}. 


(6.17) 


Equation  (6.17)  implies  that  the  operator  Sh  is  positive  definite.  Thus,  the  CFL  condition 
assures  the  stability  of  the  scheme  (4.11).  ■ 
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7  Construction  of  Perfectly  Matched  Layers 


We  use  a  perfectly  matched  layer  model  that  is  derived  in  [4],  The  construction  of  this 
model  follows  the  derivation  of  Sacks  et  al.  in  [25].  We  outline  the  salient  features  of  this 
construction  below.  We  begin  with  a  form  of  Maxwell’s  equations  suitable  for  general  media 
which  permit  both  electric  and  magnetic  currents  but  do  not  contain  unbalanced  electric 
charges,  given  as 


dB 

dt 

—VxE  —  JM  ; 

(Faraday’s  Law) 

<9D 

~dt 

VxH  -  JE  ; 

(Ampere’s  Law) 

(7.1) 

V  B  = 

0  ; 

(Gauss’s  Law  for  the  magnetic  field) 

V  D  = 

0  ; 

(Gauss’s  Law  for  the  electric  field) 

Constitutive  relations  which  relate  the  electric  and  magnetic  fluxes  (D,B)  and  the  electric 
and  magnetic  currents  (Je,  Jm)  to  the  electric  and  magnetic  fields  (E,  H)  are  added  to  these 
equations  to  make  the  system  fully  determined  and  to  describe  the  response  of  a  material 
to  the  electromagnetic  fields.  In  empty  space,  these  constitutive  relations  are  D  =  e0E,  and 
B  =  /j0H,  and  Je  =  Jm  =  0,  where  eo  and  /To  are  the  permittivity  and  the  permeability  of 
free  space.  In  general,  there  are  different  possible  forms  for  these  constitutive  relationships. 
In  a  frequency  domain  formulation  of  Maxwell’s  equations,  these  can  be  converted  to  linear 
relationships  between  the  dependent  and  independent  quantities  with  frequency  dependent 
coefficient  parameters. 

We  will  derive  a  PML  model  in  the  frequency  domain  and  then  obtain  a  PML  model  in 
the  time  domain  by  taking  the  inverse  Fourier  transforms  of  the  frequency  domain  equations. 
To  this  end,  we  consider  the  time-harmonic  form  of  Maxwell’s  equations  (7.1)  given  by 

\ujB  =  —VxE  —  J  M, 

iwD  =  VxH  —  Je, 

(7-2) 

V  •  B  =  0, 

,  V-D  =  0, 


where  for  every  field  vector  V,  V  denotes  its  Fourier  transform,  and  we  have  the  constitutive 


laws 

B  =  [/i]H, 

D  =  [e]E, 

j  M  =  [cm]H, 

J  E  =  [<7e]E. 


(7.3) 
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Here,  the  square  brackets  indicate  a  tensor  quantity.  Note  that  when  the  density  of  electric 
and  magnetic  charge  carriers  in  the  medium  is  uniform  throughout  space,  then  V  •  J#  =  0 
and  V  •  Jm  =  0.  We  define  the  new  tensors 


\tA  —  l/'l  + 


e  =  e  + 


Wm] 

lev 

We] 
i  u) 


(7.4) 


Using  the  definitions  (7.4)  we  define  two  new  constitutive  laws  that  are  equivalent  to  (7.3), 
given  by 


Br 


Dr 


=  MH, 

=  WE. 


(7.5) 


Using  (7.5)  in  (7.2),  Maxwell’s  equations,  in  time-harmonic  form,  become 

(  icuBnew  =  — VxE, 
icuDneW  =  VxH, 

V  •  Bnew  =  0, 

V  •  Dnew  =  0. 


(7.6) 


To  apply  the  perfectly  matched  layer  to  electromagnetic  computations,  we  surround  the 
computational  domain  with  layers  of  finite  depth  in  which  the  outgoing  waves  are  trapped 
and  attenuated.  These  layers  are  backed  with  a  conventional  boundary  condition,  such 
as  a  perfect  electric  conductor  (PEC).  This  truncation  of  the  layer  will  lead  to  reflections 
generated  at  the  PEC  surface,  which  can  propagate  back  through  the  layer  to  re-enter  the 
computational  region.  In  this  case,  the  reflection  coefficient  R,  is  a  function  of  the  angle  of 
incidence  6*,  the  depth  of  the  PML  6,  as  well  as  the  specific  form  of  the  tensors  [e]  and  [fi] . 
As  shown  in  [25],  the  constitutive  laws  for  the  PML  region  are 


[e]  =  e„[S],  (7.7) 

[P]  =  ho  [A],  (7.8) 

and  the  correct  form  of  the  tensor  [A]  which  appears  in  the  constitutive  laws  is  the  product 

[S]  =  (7-9) 

in  which  component  [S']a  in  the  product  in  (7.9)  is  responsible  for  attenuation  in  the  a 
direction,  for  a  —  x,y,  z.  All  three  of  the  component  tensors  in  (7.9)  are  diagonal  and  have 
the  forms 


’  s*1 

0 

0  " 

[  sy 

0 

0  " 

Sz 

0  0 

[S]*  = 

0 

0 

;  [«]»  = 

0  ^ 

0 

;  [S\,= 

0 

0 

0 

0 

L 0 

0 

Sy  _ 

_  0 

0  sj1  _ 
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The  constitutive  laws  for  the  perfectly  matched  layer  are  now  given  to  be 


(  B, 


=  //o[S]  H, 


Dnew  —  6o[aS]E. 


(7.11) 


The  parameters  sx,sy,  and  sz  in  the  PML  layers  are  chosen  in  order  for  the  attenuation  of 
waves  in  the  PML  to  be  sufficient  so  that  the  waves  striking  the  PEC  surface  are  negligible 
in  magnitude.  Perfectly  matched  layers  are  then  placed  near  each  edge  (face  in  3D)  of  the 
computational  domain  where  a  non- reflecting  condition  is  desired.  This  leads  to  overlapping 
PML  regions  in  the  corners  of  the  domain.  When  designing  PML’s  for  implementation,  it 
is  important  to  choose  the  parameters  sa  so  that  the  resulting  frequency  domain  equations 
can  be  easily  converted  back  into  the  time  domain.  The  simplest  of  these  which  we  employ 
here  [13]  is 

sa  =  1  -| - — ,  where  aa  >  0  a  =  x,y,z.  (7.12) 

icue0 

Gedney  [13]  suggests  a  conductivity  profile 

t  n,  Chnaxl®  tto|  1 

^ - ;  a  =  x,y,z.  (7.13) 

where  5  is  the  depth  of  the  layer,  a  =  a0  is  the  interface  between  the  PML  and  the  com¬ 
putational  domain,  and  m  is  the  order  of  the  polynomial  variation.  Gedney  remarks  that 
values  of  m  between  3  and  4  are  believed  to  be  optimal.  For  the  conductivity  profile  (7.13), 
the  PML  parameters  can  be  determined  for  given  values  of  m,  S,  and  the  desired  reflection 
coefficient  at  normal  incidence  R0,  as 


(m  +  1)  ln(l/i?0) 
CTmax  ~  2^5  ’ 


(7.14) 


rj  being  the  characteristic  wave  impedance  of  the  PML.  Empirical  testing  suggests  that,  for 
a  broad  range  of  problems,  an  optimal  value  of  amax  is  given  by 


•^opt  ~ 


m  +  1 

1507rAav/e7’ 


(7.15) 


where  Aa  is  the  space  increment  in  the  a  direction  and  er  is  the  relative  permittivity  of  the 
material  being  modeled.  In  the  case  of  free  space  er  =  1. 

From  the  time-harmonic  Maxwell’s  curl  equations  in  the  PML  (7.6)  and  (7.11),  Ampere’s 
and  Faraday’s  laws  can  be  written  in  the  most  general  form  as 


f  iuy0[S]K 


-VxE 

VxH 


v 


18 


(Faraday’s  Law) 
(Ampere’s  Law) 


(7.16) 


In  the  presence  of  the  diagonal  tensor  [S'],  a  plane  wave  is  purely  transmitted  into  the  PML. 

To  obtain  the  two  dimensional  model  of  the  PML,  we  assume  no  variation  in  the  z 

d 

direction  (i.e.,  —  =  0).  In  the  two  dimensional  TM  mode  the  electromagnetic  field  has 
dz 

three  components,  Ez,  Hx,  and  Hy.  In  this  case,  we  have  az  =  0  and  sz  =  1  in  the  PML  and 
the  time-harmonic  Maxwell’s  equations  (7.16),  in  the  PML  medium  can  be  written  in  scalar 
form  as: 


icu  /jjQ  Hx 

Sry. 


ic ojjLQ—Hy 


0EZ 
'  dy  ’ 
dEz 


dx  ’ 


(7.17) 


iuj€oSxSyEz 


OH,,  OH,. 


dx  dy 

To  avoid  a  computationally  intensive  implementation,  we  do  not  insert  the  expressions 
for  sx,  sy  and  sz,  obtained  via  (7.12),  into  (7.16),  and  transform  to  the  time  domain.  Instead, 
we  define  suitable  constitutive  relationships  that  facilitate  the  decoupling  of  the  frequency 
dependent  terms  [27].  To  this  end,  we  introduce  the  fields 


Bx  yo^x  Hx, 

Ey  =  yo Sy  Hy, 
Dz  =  y0syEz. 


(7.18) 


Substituting  the  definitions  (7.18)  in  (7.17),  using  the  defining  relations  for  sx  and  sy  from 
(7.12),  and  then  transforming  into  the  time  domain  by  using  the  inverse  Fourier  transform, 
yields  an  equivalent  system  of  time-domain  differential  equations,  which  is  the  two  dimen¬ 
sional  TM  mode  of  the  PML  given  as 


dB 

dt 

dH 

~dt 

dD 

dt 

dE 

<  dt 


— — £2B  - 
1  dB 
Ho  dt 

- &xD  + 

- ( 7yE  + 

e0 


curl  E, 
— SiB, 

6q/J,q 

curl  H, 

1  dD 
eo  dt  ’ 


with  H  =  (Hx,  Hy)T ,  B  =  (. Bx ,  By)T ,  E  =  Ez  and  D  =  Dz.  Also, 


Si 


0  y  y  =  (  av  0  \ 

0  Vy  )'  2  {  0  axJ- 


(7.19) 


(7.20) 


Thus,  the  PML  model  consists  in  solving  system  (7.19)  for  the  six  variables,  Bx,  By,  Hx, 

Hy,  Dz,  Ez. 
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In  order  to  use  the  PML  model,  instead  of  the  first  order  Silver-Miiller  boundary  condi¬ 
tion,  we  replace  the  discrete  model  (5.1)- (5.3)  by  the  PML  model  and  then  account  for  the 
Dirichlet  condition  on  the  boundary  of  the  obstacle  using  (5.6).  The  distributed  multiplier 
formulation  with  perfectly  matched  layers  is  then  given  to  be 

Find  (E%ee  |n+1,  Dl+\  H.nh+1*  ,Bnh+1*)  EUhxUhxVhxVh  such  that  for  all  Vh  e  Vh,  for  all 
4*h  ^ 


(i)  = 

(ii)  *h)  = 


,n+i 


e0 

ho 

1  ,  ^n+i 


-  (curl  E]';  ,  *h), 


+ 


1  (SrB 


6o/iQ 


(hi)  (DAtD,  \<j>h)  = - ',4)  +  (curl<^,Hfc 

<4) 


■n+i, 


(iv)  (DMEjr\n+i,M  =  --kc + 

Co 


e0 


-(DmD 


n+i 


In  the  above 


and 


DmBT\ 


vhi/2  _  ^yeT+1  -  K 


At 


Ke  I 


H-i/2  =  ^rr+1 
2 


(7.21) 


(7.22) 


(7.23) 


Once  we  obtain  the  solution  to  the  system  (7.21),  the  solution  to  the  scattering  problem  is 
obtained  by  solving  system  (5.6). 


8  Numerical  Examples:  Scattering  by  a  Disk 

We  consider  the  scattering  of  the  harmonic  planar  waves  e-1(p*-k 'x^  by  a  perfectly  reflecting 
disk  whose  radius  is  0.25  m.  The  disc  is  located  at  the  center  of  the  domain  [0,  3.5]  x  [0,3.5], 
The  frequency,  /,  is  0.6  GHz,  and  the  wavelength,  L,  is  0.5  m.  The  angular  frequency  is 
p  =  27 rf.  The  wave  illuminates  c o  from  the  left  and  propagates  horizontally.  The  distance 
from  the  disk  to  the  absorbing  boundary  is  3  wavelengths.  We  have  used  a  rectangular  mesh 
consisting  of  113  x  113  nodes,  with  the  mesh  step  size  h  =  0.5/16  m.  The  time  step  is 
At  =  27t/(25 p).  Thus,  the  Courant  number  is  (cAt)/h  =  0.64.  We  have  also  considered 
mesh  refinements  in  order  to  estimate  the  accuracy  of  our  solution. 

In  Figure  4  we  plot  the  number  of  degrees  of  freedom  (DOF)  of  the  Lagrange  multiplier 
on  the  boundary  of  the  disk,  du,  as  a  function  of  the  mesh  ratio,  hg^/h,  for  different  dis¬ 
cretizations,  where  hgu  is  defined  as  the  step  size  on  the  boundary  of  the  disk.  As  can  be  seen 
from  Figure  4,  for  fine  meshes,  as  opposed  to  coarse  meshes,  bigger  changes  in  the  DOF  on 
the  boundary  of  the  disk  result  from  a  small  change  in  the  mesh  ratio.  The  Uzawa  algorithm 
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Figure  4:  The  number  of  degrees  of  freedom  (DOF)  of  the  Lagrange 
multiplier  on  the  boundary  of  the  disk  versus  the  mesh  ratio  hg^/h. 


converges  in  a  finite  number  of  iterations  for  values  of  hgu/h  between  1.5  and  3.  However, 
for  certain  values  between  1.5  and  1.8  the  behavior  of  the  Uzawa  algorithm  is  unstable  with 
respect  to  number  of  iterations.  Thus  we  consider  values  for  hg^/h  between  1.8  and  3. 

For  this  test  problem  the  exact  solution  is  known  when  T  is  located  at  infinity.  In  Figures 
5,  6  and  7,  we  plot  the  error  (point-wise  difference)  between  each  computed  solution  and  the 
exact  solution  for  discretizations  with  16,  32  and  64  nodes  per  wavelength,  respectively.  It  is 
clearly  observed  that,  as  the  mesh  is  refined,  the  error  in  the  fictitious  domain  method  with 
the  Silver-Miiller  condition  is  dominated  by  reflections  from  the  artificial  boundary.  In  the 
case  of  the  PML  model  the  error  in  the  discretization  of  the  Lagrange  multiplier  dominates 
the  total  error. 

Next,  we  define  the  relative  error  (RE)  between  the  exact  solution  and  a  computed 
solution  as 

-p-p  _  I  Inexact  —  Ec\\l2(CI) 

-EexactU  l2(Q) 

where  Eex act  stands  for  the  exact  solution,  and  E(-  denotes  a  computed  solution.  In  Figure 
8  we  plot  the  relative  error  for  the  fictitious  domain  method  with  a  PML  of  thickness  L/4, 
against  the  mesh  ratio  hg^/h,  for  different  discretizations.  From  Figure  8  (left)  we  can  see 
that  the  relative  error  can  vary  by  a  factor  of  2  for  different  values  of  the  mesh  ratio.  In 
this  figure  we  also  plot  (right)  ratios  of  relative  errors  between  successive  mesh  refinements 
obtained  from  Figure  8  (left).  The  solid  line  represents  the  ratio  of  the  relative  error  for  a 
discretization  with  16  nodes  per  wavelength,  and  the  relative  error  of  a  discretization  with 


8.1) 
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Error 


Error 


Figure  5:  Plot  of  the  error  between  the  exact  solution  and  the  fictitious 
domain  method  for  a  discretization  with  16  nodes  per  wavelength:  (left) 
with  the  Silver-Miiller  boundary  condition;  (right)  with  a  4  cell  PML. 


Error 


Error 


Figure  6:  Plot  of  the  error  between  the  exact  solution  and  the  fictitious 
domain  method  for  a  discretization  with  32  nodes  per  wavelength:  (left) 
with  the  Silver-Miiller  boundary  condition;  (right)  with  a  4  cell  PML. 
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Error 


Error 


Figure  7:  Plot  of  the  error  between  the  exact  solution  and  the  fictitious 
domain  method  for  a  discretization  with  64  nodes  per  wavelength:  (left) 
with  the  Silver-Miiller  boundary  condition;  (right)  with  a  4  cell  PML. 


32  nodes  per  wavelength.  Similarly,  the  dashed  line  denotes  the  ratio  of  relative  errors  for 
discretizations  with  32  and  64  nodes  per  wavelength.  Again,  we  observe  that  the  ratios  can 
vary  by  almost  a  factor  of  3. 

In  Figures  9  we  plot  the  maximum  and  the  minimum  iteration  counts,  respectively,  that 
are  required  for  the  convergence  of  the  Uzawa  algorithm,  as  a  function  of  the  mesh  ratio. 
These  results  are  for  the  fictitious  domain  method  with  a  PML  of  thickness  L/4.  The  number 
of  iterations  for  the  three  discretizations  is  seen  to  be  bounded  by  20. 

In  Table  1  we  present  relative  errors  and  maximum  and  minimum  iteration  counts  for 
the  fictitious  domain  method  with  PML’s  of  thickness,  L/4,  L/2  and  L.  The  relative  error 
in  all  three  cases  is  almost  the  same.  Thus,  we  do  not  obtain  any  benefit  by  increasing  the 
thickness  of  the  PML  as  the  dominating  error  is  due  to  the  discretization  of  the  Lagrange 
multiplier.  This  can  also  be  observed  in  the  Figures  5,  6,  and  7.  A  quarter  wavelength  thick 
PML  is  sufficient  to  obtain  significant  improvements  over  the  first  order  boundary  condition. 

In  Table  2  we  present  relative  errors  and  maximum  and  minimum  iteration  counts  for 
the  fictitious  domain  method  with  the  Silver-Miiller  boundary  condition.  As  observed  in 
Figures  5,  6,  and  7,  reflections  from  the  artificial  boundary  dominate  the  error.  Thus,  we  do 
not  expect  to  see  much  improvement  in  the  error  as  the  mesh  is  refined.  As  seen  in  Table  1, 
the  relative  error  for  a  discretization  with  64  nodes  per  wavelength  is  4  times  smaller  than 
the  error  for  the  same  discretization  with  the  Silver-Miiller  boundary  condition. 
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Maximum  Iterations  in  Uzawa  Relative  erro 


Figure  8:  (left)  Plot  of  the  relative  error  versus  the  mesh  ratio  hg^/h  for 
three  different  discretizations,  (right)  Plot  of  ratios  of  successive  relative 
errors  versus  the  mesh  ratio  hgujh. 


Figure  9:  The  maximum  (left)  and  minimum  (right)  number  of  iterations 
required  for  the  Uzawa  algorithm  versus  the  mesh  ratio  hgu/h  for  three 
different  discretizations. 
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PML 

h 

At 

RE 

Max  Iter 

Min  Iter 

L/4 

L/ 16 

L/(25c) 

5.867e-2 

14 

12 

L/ 32 

L/(50c) 

2.389e-2 

13 

12 

L/64 

L/(100c) 

9.830e-3 

13 

12 

L/2 

L/16 

L/(25c) 

5.815e-2 

14 

12 

L/ 32 

L/(50c) 

2.389e-2 

13 

12 

L/64 

L/(100c) 

9.830e-3 

13 

12 

L 

L/16 

L/(25c) 

5.815e-2 

14 

12 

L/ 32 

L/(50c) 

2.389e-2 

13 

12 

L/64 

L/(100c) 

9.829e-3 

13 

12 

Table  1:  Table  of  relative  errors  of  the  fictitious  domain  solutions,  for 
PML’s  of  varying  thickness,  computed  with  respect  to  the  exact  solution. 


h 

At 

SM 

RE 

Max  Iter 

Min  Iter 

L/16 

L/ (25c) 

7.026e-2 

14 

12 

L/32 

L/ (50c) 

4.519e-2 

13 

12 

L/64 

L/ (100c) 

3.854e-2 

13 

12 

Table  2:  Table  of  relative  errors  of  the  fictitious  domain 
solution,  with  the  Silver-Miiller  (SM)  boundary  condition, 
computed  with  respect  to  the  exact  solution  for  different 
discretizations. 


Finally,  in  Table  3,  we  present  relative  errors  and  maximum  and  minimum  iteration 
counts  for  the  fictitious  domain  method  with  a  PML  of  thickness  L/ 4,  for  different  values 
of  the  mesh  ratio  kg^/h,  and  for  different  discretizations.  For  a  given  discretization,  the 
relative  errors  are  comparable  for  different  values  of  the  mesh  ratio.  We  observe  that  the 
ratios  between  successive  relative  errors,  for  a  fixed  value  of  the  mesh  ratio,  is  approximately 
2.  Thus,  the  spatial  accuracy  of  the  fictitious  domain  method,  based  on  these  results,  seems 
to  be  about  first  order. 
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h 

At 

hdco/h 

PML  (L/4) 

RE 

Max  Iter 

Min  Iter 

L/ 16 

LI  (25c) 

2.0 

5.867e-2 

14 

12 

2.4 

7.577e-2 

3 

3 

2.7 

5.412e-2 

5 

5 

L/ 32 

L/ (50c) 

2.0 

2.389e-2 

13 

12 

2.4 

2.312e-2 

14 

13 

2.7 

1.691e-2 

14 

13 

L/64 

LI  (100c) 

2.0 

9.830e-3 

13 

12 

2.4 

9.355e-3 

15 

13 

2.7 

6.641e-3 

14 

12 

L/128 

L/(200c) 

2.0 

6.972e-3 

16 

15 

2.4 

4.957e-3 

16 

14 

2.7 

4.679e-3 

13 

11 

Table  3:  Table  of  relative  errors  of  the  fictitious  domain  solutions  for  a 
PML  of  thickness  L/4,  computed  with  respect  to  the  exact  solution,  for 
different  values  of  the  mesh  ratio  hg^/h  and  different  discretizations. 


In  Table  4  we  compare  the  fictitious  domain  approach  to  a  staircase  approach  using  the 
finite  difference  time  domain  method.  As  can  be  seen  from  the  table,  the  fictitious  domain 
method  provides  a  significant  improvement  over  the  staircase  approximation.  This  is  also 
evident  from  Figures  10  and  11,  which  compare  the  errors  for  both  methods  for  16  and  64 
nodes  per  wavelength. 


N 

L/h 

Staircase 

Fictitious  Domain 

U32 

16 

1.959e-l 

5.867e-2 

2252 

32 

9.997e-2 

2.389e-2 

4492 

64 

4.871e-2 

9.830e-3 

8972 

128 

2.619e-2 

6.972e-3 

Table  4:  Table  of  relative  errors  for  the  fictitious  domain  solution 
for  a  PML  of  thickness  L/4,  and  relative  errors  for  a  staircase 
approximation  for  different  nodes  per  wavelength. 


So  far  we  have  presented  results  in  which  the  frequency  /,  and  hence  the  wavelength  L, 
in  the  domain  was  fixed.  As  the  frequency  is  increased,  i.e.,  the  wavelength  is  decreased, 
the  effects  of  dispersion  start  to  degrade  the  solution.  The  error  in  the  solution  is  no  longer 
dominated  by  the  error  in  the  discretization  of  the  Lagrange  multiplier.  The  error  at  higher 
frequencies  is  dominated  by  large  phase  errors,  which  accumulate  over  time  and  can  sig¬ 
nificantly  affect  the  solution.  To  study  the  errors  that  arise  at  high  frequencies,  we  have 
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calculated  the  relative  errors  in  the  case  when  /  =  1.2, 2.4  and  4.8  GHz.  The  relative  error  is 
a  combination  of  the  error  in  the  amplitude  of  the  solution  as  well  as  the  error  in  the  phase. 
To  see  the  dominance  of  the  phase  error,  we  also  calculate  a  relative  amplitude  error  and  a 
phase  error  for  each  frequency.  In  these  calculations  the  size  of  the  domain  is  fixed.  We  use 
a  fixed  step  size  h  =  0.5/128,  which  uses  128  nodes  per  wavelength  at  the  lowest  frequency 
of  0.6  GHz,  on  the  square  domain  [0,3.5]  x  [0,3.5].  Thus,  we  have  897  x  897  nodes  with 
N  =  897.  The  time  step  is  chosen  so  that  the  stability  condition  as  before  is  t]  =  0.64.  The 
results  are  all  computed  after  1400  iterations.  All  results  are  computed  using  a  4  cell  PML 
and  the  fictitious  domain  method.  We  do  not  observe  any  significant  reduction  in  the  errors 
by  increasing  the  thickness  of  the  PML  layer. 

In  Figure  12  we  plot  a  top  view  of  the  solutions  with  /  =  0.6  GHz  (top)  and  /  =  1.2 
GHz  (bottom).  In  Figure  13  we  plot  a  top  view  of  the  solutions  with  /  =  2.4  GHz  (top)  and 
/  =  4.8  GHz  (bottom).  The  computed  solutions  qualitatively  compare  well  with  the  exact 
solution. 

Let  Rex act  and  /exact  denote  the  real  and  imaginary  parts  of  the  exact  solution.  Similarly, 
let  Rc  and  Ic  denote  the  real  and  imaginary  parts  of  our  computed  solution.  We  define  the 
phase  error  (PE)  as 


PE(x)  =  tan  1 


f  Axact  (-X)  \ 
\  Rex  act  (-X)  / 


—  tan  1 


(8.2) 


We  also  calculate  the  phase  error  in  degrees  per  node  as 


Phase  error  = 


360 
27 tN 


N 2 

ElpEMI 


,fc=l 


1/2 

degrees/node, 


(8.3) 


where  Xk  is  a  node  of  the  finite  element  triangulation  7^.  The  amplitude  error  (AE)  is  defined 
as 

AE(X)  =  y/(Iexact (Z))2  +  (i?exact(x))2  -  y/(/C(x))2  +  (Rc(z))2-  (8.4) 

Next,  we  calculate  a  relative  amplitude  error  (RAE)  for  each  frequency  which  will  be  the 
ratio  of  the  L2  norm  of  the  amplitude  error  to  the  L2  norm  of  the  amplitude  of  the  exact 
solution  in  the  computational  domain.  In  Figure  14  we  plot  linear  gray  scale  images  of  the 
phase  error  (top)  and  the  amplitude  error  (bottom)  over  the  square  domain.  In  this  figure 
we  can  see  that  the  phase  error  is  the  smallest  along  the  grid  diagonals  and  it  is  the  largest 
along  the  axis  of  the  mesh. 

In  Table  5  we  present  the  (total)  relative  errors  for  the  real  and  imaginary  parts  of  the 
solution.  As  expected  the  relative  errors  increase  as  the  frequency  is  increased.  We  also 
compare  the  relative  amplitude  error  RAE  which  is  calculated  using  the  amplitude  error  AE 
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Fictitious  Domain  Approximation 


Staircase  Approximation 


Figure  10:  (left)  Plot  of  the  error  between  the  exact  solution  and  the 
fictitious  domain  method  with  a  4  cell  PML.  (right)  Plot  of  the  error 
between  the  exact  solution  and  a  staircase  approximation.  In  both  cases 
we  use  a  discretization  with  16  nodes  per  wavelength. 


Fictitious  Domain  Approximation 


Staircase  Approximation 


Figure  11:  (left)  Plot  of  the  error  between  the  exact  solution  and  the 
fictitious  domain  method  with  a  4  cell  PML.  (right)  Plot  of  the  error 
between  the  exact  solution  and  a  staircase  approximation.  In  both  cases 
we  use  a  discretization  with  64  nodes  per  wavelength. 
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3.5 


Frequency:  0.6  GHz 


Frequency:  1.2  GHz 


3.5 


0.5  0.5 

o- . -  o- . 

0  0.5  1  1.5  2  2.5  3  3.5  0  0.5  1  1.5  2  2.5  3  3.5 


x  axis  x  axis 

Figure  12:  A  top  view  of  the  computed  solution  (real  part)  for  a  harmonic 
planar  wave  with  frequency  /  =  0.6  GHz  and  L  =  0.5m  (left),  and  for 
a  harmonic  planar  wave  with  frequency  /  =  1.2  GHz  and  L  =  0.25m 
(right). 


defined  in  (8.4).  Finnally  we  compare  the  phase  error  in  degrees  per  node,  defined  in  (8.3), 
for  each  case.  We  observe  that  the  relative  amplitude  error  is  significantly  better  than  the 
total  relative  errors.  As  can  be  seen  from  the  table  the  phase  error  increases  significantly  at 
higher  frequencies.  For  /  =  0.6  GHz,  the  phase  error  is  0.37  degrees  per  node.  This  error 
increases  to  15.69  degrees  per  node  when  /  =  4.8  GHz. 


/  (GHz) 

L  (m) 

L/h 

RE  (Real) 

RE  (Imag) 

RAE 

Phase 

0.6 

0.5 

128 

6.97e-3 

7.77e-3 

3.31e-3 

0.37 

1.2 

0.25 

64 

1.57e-2 

1.52e-2 

5.72e-3 

0.73 

2.4 

0.125 

32 

4.97e-2 

4.72e-2 

l.lle- 2 

2.29 

4.8 

0.0625 

16 

2.89e-l 

2.93e-l 

2.57e-2 

15.69 

Table  5:  Table  of  errors  for  the  fictitious  domain  solution  for  a  PML  of 
thickness  L/4,  at  different  frequencies.  The  relative  error  for  the  real  and 
imaginary  parts  of  the  solution  is  given.  RAE  is  a  relative  amplitude  error 
and  the  phase  error  in  degrees  per  node  in  each  case  is  provided. 
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y  axis 


Frequency:  2.4  GHz  Frequency:  4.8  GHz 
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x  axis  x  axis 


Figure  13:  A  top  view  of  the  computed  solution  (real  part)  for  a  harmonic 
planar  wave  with  frequency  /  =  2.4  GHz  and  L  =  0.125m  (left),  and  for 
a  harmonic  planar  wave  with  frequency  /  =  4.8  GHz  and  L  =  0.0625m 
(right). 
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Figure  14:  A  linear  gray  scale  image  of  the  phase  error  in  radians  (left) 
and  the  amplitude  error  (right)  over  the  square  domain  [0,  3.5]  x  [0,3.5], 
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9  Conclusions  and  Future  Work 


In  this  paper  we  a  have  applied  a  distributed  Lagrange  multiplier  based  fictitious  domain 
method  to  the  solution  of  the  two  dimensional  Maxwell’s  equations  in  the  exterior  of  a 
domain  u>  with  a  perfectly  conducting  condition  on  the  boundary  of  oj.  Such  a  method  was 
applied  to  the  solution  of  the  two  dimensional  wave  equation  in  [5].  Thus,  in  this  paper  we 
have  extended  the  application  of  the  distributed  multiplier  to  electromagnetic  waves.  The 
idea  behind  the  fictitious  domain  method  is  to  extend  the  electromagnetic  solution  inside  the 
obstacle  and  enforce  the  perfectly  conducting  condition  on  the  boundary  of  the  obstacle  via 
the  introduction  of  a  distributed  Lagrange  multiplier.  This  distributed  multiplier  is  defined 
on  the  boundary  of  the  obstacle  as  well  as  in  the  interior  of  the  obstacle. 

After  presenting  our  fictitious  domain  formulation,  we  have  derived  energy  identities  to 
demonstrate  the  wellposedness  and  stability  requirements  of  the  method.  An  interesting  fact 
is  that  the  stability  condition  is  the  same  as  in  the  case  of  the  problem  without  an  obstacle. 
We  have  presented  numerical  computations  which  show  that  our  fictitious  domain  method 
is  more  accurate  than  the  finite  difference  approach.  Even  though  the  method  remains  first 
order  accurate  with  respect  to  the  mesh  step  size,  the  error  is  much  better  as  compared  to 
the  staircase  approximation  of  a  the  corresponding  finite  difference  scheme.  Thus,  this  paper 
provides  a  simple  and  much  more  accurate  alternative  to  the  popular  FDTD  method. 

We  have  focused  our  study  on  the  two  dimensional  TM  mode  of  Maxwell’s  equations.  In 
the  future  we  will  extend  this  approach  to  the  TE  mode  in  two  dimensions  as  well  as  to  the 
full  three  dimensional  Maxwell’s  equations. 
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